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ABSTRACT 

In this paper we study uniform high order spectral methods to solve multi-dimensional 
Euler equations for gas dynamics. Uniform high order spectral approximations with spectral 
accuracy in smooth regions of solutions are constructed by introducing the idea of the Essen- 
tially Non-Oscillatory (ENO) polynomial interpolations into the spectral methods. Based 
on the new approximations, we propose nonoscillatory spectral methods which possess the 
properties of both upwinding difference schemes and spectral methods. We present numerical 
results for the inviscid Burgers’ equation, and for one dimensional Euler equations including 
the interactions between a shock wave and density disturbance, Sod’s and Lax’s shock tube 
problems, and the blast wave problem. Finally, we simulate the interaction between a Mach 
3 two dimensional shock wave and a rotating vortex. 
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1 Introduction 


Recently, high order numerical methods have attracted considerable attention^ for the sim- 
ulation of flows with shock waves and different scales, especially for turbulent flows affected 
by shock wave interactions. Those high order methods are expected to produce nonoscil- 
latory sharp shock profiles without excessive overall numerical diffusion and, at the same 
time, be able to resolve the small scales of the flow field elsewhere. Recent results with 
essentially nonoscillatory (ENO) finite difference methods have shown considerable progress 
in this direction [9], [17]. Spectral methods, as high order global methods, have been very 
successful in the study of turbulent flows and flow transition problems when the solutions of 
the fluid problems are smooth. For those problems, spectral methods have been shown to 
have an accuracy higher than any algebraic order (so called spectral accuracy) [5]. However 
it remains to show that spectral methods will also be successful in computing flows with 
shock waves. 

In this paper, we continue our previous work [4] in designing essentially nonoscillatory 
spectral methods for computing the weak solutions of the hyperbolic system of conservation 
laws 


u t + f(u)« + g(u) v = 0 (1.1) 

u(x,y,0) = u 0 (a;,y). (1.2) 

Here, as usual, u = (tii, • ■ • ,u,) T is a state vector and f(u),g(u) are the vector- valued 
flux functions of s components. The system is assumed to be hyperbolic in the sense that 
for any real vector £ = (£i, £ 2 ), the matrix always has s - real eigenvalues and 

a complete set of eigenvectors. The solutions to (1.1) usually develop discontinuities in the 
form of shock waves and contact discontinuities. 

In applying spectral methods to problems having discontinuous solutions, a key issue is 
how to deal with the Gibbs phenomenon caused by the discontinuities of the solutions. The 
overall accuracy of spectral methods will be, at most, first order everywhere in the presence 
of Gibbs oscillations. There are various filtering techniques to recover spectral accuracy in 
regions away from the discontinuities [8], [14]. On the other hand, one-sided filtering can 
be used to obtain uniform convergence in regions close to the discontinuities [2]. As another 
approach to treat the Gibbs oscillation, in [4] we proposed a nonoscillatory spectral approxi- 
mation to discontinuous solutions by adding piecewise linear functions, such as sawtooth-like 
functions and step functions, to the conventional Fourier trigonometric or Chebyshev poly- 
nomial spaces. Those additional functions are used to resolve the discontinuities in the 
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solutions caused by shock waves and contact discontinuities. The cell- averaged form of 
(1.1) is used to formulate the numerical schemes, resulting in Godunov-type shock capturing 
algorithms. The usual reconstruction step between cell averages and point values of the nu- 
merical solutions in such schemes can be performed efficiently with Fast Fourier Transforms. 
However, a common problem with cell-averaged formulation is the costly implementation of 
the reconstruction in multi-dimensional problems. 

In this paper, we adopt the same philosophy as in [4], however, a more robust and 
sophisticated technique will be introduced. With the new technique, we will be able to 
achieve global convergence up to any given m-th order (m > 0) and , meanwhile, retain 
spectral accuracy in the regions away from the discontinuities. In order to achieve these goals, 
we incorporate the main idea of the ENO polynomial interpolations [9] into our construction 
of uniform spectral approximations. We also introduce the idea of upwind differencing from 
conservative finite difference methods into the design of the spectral schemes. The idea 
of upwinding has proven very successful in capturing shock wave fronts and producing the 
right entropy solutions. By using local Riemann solvers and flux limiters, modern shock 
capturing finite difference schemes, such as the TVD schemes [9] , MUSCL type schemes 
[19], FCT schemes [1], and the more recent ENO schemes [9] [17], produce very satisfactory 
shock profiles and entropy satisfying solutions. The nonoscillatory spectral approximations 
proposed in this paper will enable us to bring the upwinding idea into the framework of 
spectral methods. Meanwhile, the spectral schemes will be based directly on the conservation 
laws (l.l), not on its cell-averaged form. Thus, generalization to multi-dimensional cases is 
straightforward. 

For systems of conservation laws, in order to achieve sharp shock profile without spurious 
oscillations, numerical flux operators for the scalar equations are usually applied to the 
locally defined characteristic variables. Because of this complication, it has been realized 
that the cost of upwinding schemes is much greater than that of the centered difference 
schemes. Several attempts have been made to eliminate this shortcoming by combining 
centered difference schemes and upwinding schemes. In [13], a mixed method of centered 
difference schemes and ENO schemes was studied and, in [6], the authors suggested a type 
of nonlinear filtering technique to modify the results of the Lax-Wendroff scheme at each 
time step to produce nonoscillatory TVD solutions. Even though it is achieved in a quite 
different way from those in [6] and [13], the result in this paper will provide another example 
of blending the nice properties of both upwinding schemes and centered difference schemes 
(in this case, spectral schemes). 

This paper is organized as follows: in Section 2, we first briefly review the method 
proposed in [4], then present the new method of constructing uniform convergent, up to any 
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given ra-th order(m > 0), spectral approximations to discontinuous functions. In Section 

3, we study the nonoscillatory spectral methods for scalar conservation laws. Extensions to 
the system of conservation laws and multi- dimensional problems will be discussed in Section 

4. In Section 5, we present numerical experiments for the new methods. First, the uniform 
convergence and the spectral accuracy of the proposed spectral approximations are tested 
on discontinuous functions. Then we study the global accuracy of the spectral schemes on 
a scalar in viscid Burgers’ equation and on one dimensional Euler equations which model 
the interaction of a pure shock wave with density waves. Also we apply the scheme to 
the standard Sod’s and Lax’s test problems [16] in order to check the convergence of our 
spectral schemes to the correct entropy solutions. High order numerical results will also be 
presented for solving the interaction between two blast waves [20]. Finally we apply the 
spectral schemes to simulate the interactions between a Mach 3 two-dimensional shock wave 
and a rotating vortex. 

2 Uniform High Order Spectral Approximations 

The conventional Fourier spectral space has basis functions {e* fcx }|*|<jv- The Fourier expan- 
sions for discontinuous functions converge very slowly. For instance, consider a sawtooth-like 
function 


F'(x, x 3) 



—x 

2w — x 


for x < x l} 
for x > x„ 


( 2 . 1 ) 


where x t is the location of the discontinuity and A = ) — [ji 1 ]^ j s the jump of 

F(x,x m , A) across x,. The partial sum of the Fourier expansion of F(x,x Si A) is 


F N (x,x t) A)= £ f k (x.,A)e ik *, (2.2) 

\k\<N 

where 

fk( x s, A) = - -J o F(x, x„ A)e-' k * dx = A j ^ * L“o!’ (2.3) 

From (2.3) we can see that the Fourier coefficients f k (x t ,A) only decay like 0(1) as 
k — » oo. As a result, the convergence of (2.2) will be only first order, and moreover, the 
Gibbs oscillations near x, will be in the order of 0(1). In order to get rid of the Gibbs 
oscillations, in [4] we proposed a technique to construct essentially nonoscillatory spectral 
approximations, which we review below. 

Let u(x) be a piecewise C°° periodic function with a jump discontinuity at x, with jump 
[u] Xl and let Uff(x) be its finite Fourier expansion, the nonoscillatory spectral approximation 
is defined by 
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(2.4) 


«»(») = E «*«“* + E 

where y is an approximation of x 3 and A 1 is an approximation of [v] Xa and 

1 r 2ir 

a k = — / ti(cc)e- lAx dx. 

Z7T JO 

Since the second sum in (2.4) is actually F(x,y,A') — FN(x,y, A'), we have 

«S,(») = E K - + F(x, y, A!). (2.5) 

\k\<N 

Therefore u* N (x ) defines an approximation in the spectral space {e‘ fcx }|fc|<Ar augmented 
by sawtooth-like functions F(x,y , A'). 

The approximation defined in (2.5) yields nonoscillatory numerical results for discon- 
tinuous functions, and spectral schemes using this approximation have given high order 
accuracy for one-dimensional Euler gas dynamics equations ([2], [3]). In order for (2.5) to 
be nonoscillatory, the approximation for the location and the magnitude of the shock should 
be reasonably accurate. Second order accuracy in the location and first order accuracy in 
the magnitude are needed to ensure the uniform nonoscillatory convergence. 

In what follows, we present a different method which will be uniformly convergent up to 
any given order m > 0 and , at the same time, retain the spectral accuracy in the smooth 
regions away from the discontinuities. Furthermore, the requirement of accuracy in shock 
locations will be much relaxed and an approximation to the magnitude of the shock is not 
needed. This makes the computation more robust. 

Before we discuss the new approximation method we review two techniques to be used 
in our construction. The first one is the essentially non-oscillatory (ENO) polynomial inter- 
polation, and the second is the filtering technique for Fourier approximations. 


ENO Polynomial Interpolation 

We will follow the notation used in [9]. Let u(x) be a function defined on I = [0,27r] 
and {xj^o be the uniform mesh on 7, x» = ih,h = For simplicity of illustration, we 
assume that u(x) has only one discontinuity at x, £ I. Now given u (xj), 0 i < N , define a 
piecewise m-th order polynomial interpolant Q m (x; u) for u(x) at mesh points Xj, 0 < i < N 
as follows: 

Q m (xi]u) = u(xi) for 0 < i < N, (2.6) 

and 

Q m (x;u) = g mj . + i(x;u) for x i < x < x J+1 , (2.7) 
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where g m j+ i (x; u) is a polynomial of degree m defined below. 

Polynomial 9 m j+ i(x;u) interpolates u(x ) at (m + 1) successive points Xj, i m {j) < i < 
*m(j)+ m * The stencil of these (m+1) mesh points will be chosen according to the smoothness 
of the data u(xj) around Xj. A recursive algorithm to define i m {j) starts by defining 


= j (2.8) 

i.e. g Xj+ i(x) will be the first degree polynomial which interpolates u(x) at Xj,x J+1 . If we 
assume g fc j+ ±(x) is the k - th degree polynomial which interpolates u(x) at 

’ ' ' j (2-9) 


then we need one additional mesh point in order to define g fc+1 j + i(x). That point may be 
the nearest one to the left of stencil of (2.9) (i.e. ®t*(j)-i) or the nearest one to the right of 
the stencil of (2.9) (i.e. Xi^+fc+i)- The choice will be based on the absolute values of the 
corresponding (k + l)-th order divided differences, namely 


_ / ik(j) 1 if ■ ■ • > ®t*(j)+fc] < u [ x *n(i)> ' ’ ’ > x *k(i)+fe+i]> 

| otherwise. 


( 2 . 10 ) 


The piecewise polynomial Q m (x; u) defined in (2.6), (2.7) will give uniform nonoscillatory 
approximations to u(x) up to the discontinuities. In fact it can be shown that 


■j~^Q m {x\u) = -j-^u(x) + 0(h m+1 k ) for 0 < A; < m. (2.11) 

except for the cell containing x 4 . Using sub-cell resolutions [10] this is also correct in the 
shocked cell. 


Filtering Techniques for Fourier Approximations 

When a function u(x) is discontinuous, its Fourier approximation u# (x) will be at most 
first order everywhere [7] . However, there are several ways to recover the loss of the spectral 
accuracy in the smooth regions of the function u(x) ([8], [14]). The most common way is 
to multiply the Fourier coefficients of u(x) by a decreasing scalar factor < 7 *,. [tx* | — + 0 as 
|fc| — ► N. Then the resulting series (the filtered approximation) will be denoted by u^(x), 


u Ar( a: ) = E <r k a k e tkx . 

\k\<N 


(2.12) 


It was proven in [20] that, if cr k is derived from a scalar function cr(t), 0 < t < 1 and 
cr k = cr(^), | A: j < N , and a(t) satisfies the following conditions: 
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(2.13) 


<r(0) = 1, 
a(l) = 0, 

oW(0) = aW(l) = 0 for l<k<K, 

then Ujf(x) will converge to u{x) in the smooth regions of u(x) in the order of ). 

In the actual computations, <7* is chosen to decay exponentially in terms of the frequency 
number, 


crj t = e a ^Tr) u for |A| < N, (2.14) 

where the constant a is chosen so that o # is the machine zero and 21 is called the order of 
the exponential filtering. 

Uniform Spectral Approximations 

Now we present our new uniform high order nonoscillatory spectral approximations to 
discontinuous functions. Again, for simplicity, we assume u{x) is a periodic piecewise <7°° 
function on [0,27 t] with only one discontinuity at x,. Also we assume that the discontinuity 
has been detected within an interval [x*,xj]. 

Let us denote all the mesh points inside interval [x*, x r t ] as x t|> • • • , x^ r . Then we define a 
piecewise m - th order polynomial </?(x) which interpolates the function u(x) at mesh points 

U<i< i T , 


' q m ,j+±( x ) if * e [x j} x j+ 1] n [x l t ,x r t ] for some j, 

(p{ x) = < P|(x) if x € [0, x*], (2.15) 

. P r {x) if x E [x r „2n], 

where g mj+ i(x) was defined in (2.8) - (2.10) and Pi(x) and P P (x) are both m'-th order 
polynomials on the interval [0, x[] and [x',27 r] respectively , m' = 2m + 1, and satisfy the 
following conditions, 

Pi k) U) = 

for 0 < k < m (2.16) 

Jf’oo = <£l + j(*:). 

and 


p} k) (x r , - 2t) = 


P, ( ‘V. + 2») = 


for 0 < k < m 


(2.17) 
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Conditions (2.16) (2.17) ensure that <^(x) will be at least globally C m continuous. There 
are exact 2m + 2 = m' + 1 constraints on the m'-th polynomials Pi(x) and P r (x) respectively. 
Therefore they are uniquely defined. By (2.11), the function <p(x) will have the following 
property, 

ip(xi) = u (xi) for ij < i < i r , (2.18) 

<p(x) — u(x) = 0(h m+1 ) for x E [x*,x']. (2.19) 

Next we consider the difference between u(x) and y?(x), v(x ) = u(x) — y?(x). v(x) will 
be a C m function everywhere in [0,27r] except at x, where [u(x)] Xj = 0(h m+1 ). Moreover, 
u(xi) = 0, ii < i < i r . Therefore, the filtered Fourier interpolant Jjyv(x) will converge to 
v(x) rapidly, 


3 1 

v n( x ) = J N V ( X ) = !C a kVke xhx , 

*=-¥ 

i>k = 4 53 («(*<) “ <p(xi))e~ tkxi 
VV 1=0 

= V EM*0 - + TF EM 1 *) - 

■ iV »<»( ■ /V *>«r 

and afc is the filter in (2.14). 

Finally we define the uniform spectral approximation Vu{x) of u(x) by 
Vu{x) = <p(x) + v^(x), for x G [0,2ir]. 


Then the derivatives of u(x) will be approximated by those of Pu(x), i.e. 


d k 

dx kU 


(*) 




<i* . . 


for A; > 0. 


( 2 . 20 ) 


( 2 . 21 ) 


To see the accuracy of (2.20) to u(x), let p(x) E Co°(x',x^) be a mollifier function such 
that 


'tT 

II 

i— * 

for x near x t , 

(2.22) 

and 

if x € [x l a1 x r M ], 
otherwise. 


*•« = { :U V{X) 

(2.23) 

One can easily see that 



V*(Xi) = v(Xj) 

for 0 < i < N, 

(2.24) 
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and hence 


I N v*(x) = /jyv(®). (2.25) 

As v*(x ) € C m ( 0,2ir) and is periodic, by standard estimate it can be shown [2] that 

.IK ( ”°lk 


»*(*) - Jw”'(*)llfa 5 c ‘ 


tfm+1 > 


where c is a constant independent of N. 
On the other hand, 


by (2.19) we have 


|v*(x) - u(z)|U, = -p{x)) a v t (x)dx 

< ^£V(*)<I* 

< //>) - v ’ C *)) 2 dx > 

\\v'(x)-v(x)\\ La =0(h m+1 ). 


(2.26) 


(2.27) 


(2.28) 


It follows from (2.26) and (2.28) that 

IK*) - flW*)llt, = IK*) - ^«*(*)lk 
< |K*) - »"(*)IU, + IK(*) - J>’(*)lk (2.29) 

= 0(/i m+I ). 


Thus 


||Pu(*) - “Mlk = HW*) + fJK*)) - M*) + »(*)]||l> 

= |K*) - W*)lk = 0(/>”* +1 ), (2.30) 

which establishes the uniform (m+ 1 )- th order convergence of the approximation 'Pu of u in 
the L 2 norm. Error estimates in a higher order Sobolev norm can be derived similarly. 

The spectral convergence of Vu(x) to u(®) in the regions outside [ x \ , x T t ] follows from the 
spectral convergence of vjf(x) to u(x) in the smooth regions of v(s). 


3 Uniform High Order Spectral Methods 


In this section we study uniform high order spectral methods for conservations laws (1.1). 
First, we will consider the scalar one-dimensional conservation laws. Extensions to the 


8 



system of conservation laws and to multi- dimensional problems will be discussed in Section 


4. 

We will derive the spectral schemes using the method of lines. The time derivative and 
spatial derivatives will be discretized separately. For simplicity we only present the Euler- 
forward difference method for the time derivative. In the numerical experiments of Section 5 
a high order TVD type Runge-Kutta time discretization [17] is used. The numerical scheme 
will be written in the following conservative form, 


uf = - A(/ y+i - / M ), (3.1) 

where it" ~ u(xj,t n ),Xj = jAx,t n = nAt, and Ax and At are the spatial mesh size and the 
time step respectively, A = f j+ 1 are the numerical fluxes . 

It is observed [17] that if there is a function h(x) such that 


then 


^(“( x )) = -k L’-4 h(()d( ' 


r , , „ - h ( X >-0 

/,(«(*,)) = *-• 


(32) 


(33) 


This suggests that the numerical flux should approximate /i(x J+ i ) as Ax — ► 0. 

We construct /t(x) in the same manner as in [17] via its primitive function H(x ) modulo 
a linear function, 

H(x)= f (h(()-c)<K, (3.4) 


where c is a constant chosen so that i7(x) will be a periodic function: 

E 

j=0 


c = I K0<K = A*£/y 

J — r j=<3 

Assuming that (3.2) holds, then 

= e r kti h(od(-c(j+i)A 

L_rt X, I 


(3.5) 


(3.6) 


k=0 Jx k-$ 

= Ax ^2 f{ u ( x k )) — c (j + l)Ax for 0<j<N. 


k = 0 


We then form the uniform spectral approximation operator VH to H(x), 

VH = <p(x) + v° N (x), 


(3.7) 
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where <p(x) is the piecewise m-th polynomial defined in (2.15) and Ujy(x) is the filtered 
Fourier interpolant of data H(xj + i) — <p(xj + i ), 0 < j < N. As before, it is assumed that 
the shock discontinuity or contact discontinuity x , (for simplicity of illustration, only one 
such discontinuity is assumed to exist) has been detected in an interval [®i,a^], i.e. 



(3.8) 

If x 6 [*i,*i+i] n [x'„ X r t ] for some j, (p(x) = q m>j+ i(x). q m>j+ i(x) is 
polynomial and 

a m-th order 

9m,j+i (*<+J ) = H ( X i+$) for *m0') < i < im(j) + ™. 

(3.9) 

The stencil + m is defined recursively as in (2.8) - (2.10), however the first 

point of the stencil ii(j) is chosen according to the local Roe-speed o J+ i , 

f( u j+ 1) - /(«;) 

(3.10) 

if > 0, 

HU) -\j + 1, if a i+ x <0. 

(311) 


Then we have the following spectral algorithm 
Algorithm I (Spectral ENO-Roe) 


• Step 1, define if(x J+ 1),0 < j < N by (3.6) and their uniform spectral approximation 
VH{x) by (3.7); 

• Step 2, let 

/;+* = -^ VH ( x i+0 + c = + ^r( x ;+*) + c > ( 3 - 12 ) 

where constant c is defined in (3.5). 


Remarks 

1. The evaluation of 0 < j < N can be performed efficiently via the 

Fast Fourier Transforms with the total number of operations of order 0(N log N)] 

2. Regarding the first term, i)> 0 < j < N, only for those x J+ i € [x* , x T t ] } do 

we need to do the ENO interpolation and the numerical differentiation. If x J+ i 
is located outside [x*,xj], the derivatives can be computed analytically without 
additional computer cost; 
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3. The formal spatial accuracy of Algorithm I will be spectrally accurate in the 
smooth regions of the solution and uniformly m-order elsewhere; 

4. (Entropy Correction) the fluxes / J+ i defined in (3.12) are based on the Roe flux 
which admits “expansion shocks”. We use an entropy correction discussed in 
[17], see also Harten [ 11 ] for the origination of such corrections. 


4 Euler Equations of Gas Dynamics 


In this section we extend the scalar Algorithm I from the previous section to the system of 
Euler equations for gas dynamics for polytropic gas. With all variables in bold face denoting 
vectors, we have the following Euler equations: 


+ f(u), = 0, 
u = (p,m, E) r , 
f(u) = qu + (0 } P,qP) T , 

p = (7 - 1 )(£ - \pq 2 )- 


(4.1) 

(4.2) 

(4.3) 

(4.4) 


Here p, q , P and E are the density, velocity, pressure and total energy respectively, m — pq 
is the momentum and 7 = 1.4 is the ratio of specific heats of a polytropic gas. 


The eigenvalues of the Jacobian matrix A(u) = are 

A a (u) = q — c, A 2 (u) = q, A 3 (u) = q + c, 

where c = ( 7 _P/p )5 is the sound speed. 

The corresponding right-eigenvectors are 

/ 1 \ 


(4.5) 


ri( u ) = 


q-c 

V H-qc) 


,r 2 (u) = 



where 


»r 3 (u) = 


1 _2 


1 

q + c 
H + qc ) 


(4.6) 


H = (E + P)/p = c 1 /(- l -l)+-q\ 


is the enthalpy. 

The corresponding left eigenvectors {l*(u)} which are bi-orthonormal to {rfc(u)} in (4.6) 
are 


M u ) = ^(&2 + 9 / c >- M - 1 / c >& i )> 

1 2 ( U ) = (1 - 

H u ) = 2 ^ -9/ c >- & i9 + 1 /c, fe i)i 


(4.7) 
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where 


(4.8) 

(4.9) 


b: i = (7-1 )/c 3 , 

= \q 2 b i. 

As in the case of scalar conservation laws, we first define a vector counterpart H(x J+ i), 
0 < j < N of (3.6). The scalar quantities /(u(x*)) in (3.6) are replaced by the vectors 
f(u(x*)). We then generalize the construction of the uniform spectral approximation oper- 
ator PH to vector-valued functions in the following way. As before, the assumption about 
the shock location (3.8) still holds here. We have the uniform spectral approximation 

PH(x) = $(x) + v£(z), (4.10) 


where the components of the vector- valued function 3*(x) will be piecewise m-th order poly- 
nomials and vjy(x) will be the filtered Fourier interpolant of vector quantities H(x J+ i) — 
$(x i+ i),0 <j<N. 

»(*) will be defined separately according to whether x is inside or outside the interval 
[x*,*;]. For x 6 [xjjXj+i] n [x*,x'] for some j, first let g^ + i(s) be the ENO polynomial 
interpolant of the characteristic variables H^(x i+ t ),j — m<i<j + m. Here as usual 
the characteristic variables H^^(x i+ i ) are the projections of H(x i+ ^) on the locally defined 
characteristic fields. In this paper, we define the local characteristic fields with respect to 
the Roe-averaged state u J+ i between the states Uj,u j+ i, i.e. 

H<‘>(x, + .) = U(u i+s )H(s j+ ,) (4.11) 

for j — m <i < j + m, l < k < 3. See [15] for the definition of u J+ i . 

We thus have 


= for *»Cj) £ > £ •»(}') + m - 


(4.12) 


and, in the recursive process of choosing the stencil the first point *i(j) is determined 

according to the sign of the local eigenvalue: 

if A^(u j+ i) > 0, 
if AW(u J+ i) < 0. 


J 

i + i 


iiU) = 

We then define the vector- valued function 3>(x) as 


(4.13) 


*(®) = E £i+j(*) r *( u i+|) for x e x 3+i]- 


k— 1 


(4.14) 
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On the other hand, when x is outside of the interval [x[, x T t ], we define $(x) in the same 
way as in (2.16) and (2.17). Therefore, in those regions the components of $(x) will be 
m' — 2m + 1-th order polynomials. Globally, the components of $(x) will be C m function 
for x 6 [0, 27 r] and 

$(x i+ i) = H(x i+ i ), if x t+ i € (4.15) 

Now we present our algorithm for (4.1) 
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Algorithm II (Spectral-ENO-Roe) 


• Step 1, define vector quantities H(xy + i), 0 < j < N by (3.6) and their uniform spectral 
approximation 'PH(x) by (4.10); 

• Step 2, let 

+ c = + + c ' (416) 
where c is defined as in (3.5) with f, replaced by fy. 

Remarks 

1. All of the remarks following Algorithm I of the previous section apply to Algo- 
rithm II. However, the entropy correction is used only on the genuine nonlinear 
fields. The total number of operations in the computations of all f J+ i , 0 < j < N 
will be of order 0(N log N). Moreover, the characteristic decompositions are 
only needed for those x J+ i in [x[, x T t ] where a shock discontinuity or contact 
discontinuity has been detected; 

2. (Generalization to Mulit- dimensional Cases) For the system of conservation laws 
in two dimensions (1.1), we apply Algorithm II to f(u) and g(u) separately. 
Characteristic decompositions will be performed on their corresponding charac- 
teristic fields when needed. The same idea can be applied to higher dimensional 
cases. 


5 Numerical Results 

In this section, we will carry out several numerical experiments with Algorithm I & II. In 
implementing these algorithms, we have to choose the order of ENO interpolations in (3.7) 
or (4.10) and the interval [x\,x T ,] as defined in (3.8), which is detected to contain a shock 
or contact discontinuity. In most of the tests, we choose third order ENO interpolation, i.e. 
m = 3, unless it is mentioned otherwise. The interval [x [ , x \ ] is usually chosen to contain 6 
- 8 mesh points around a discontinuity. The numerical results show insensitivity to the size 
of [x 1 ,, x^] as long as it contains all the transition points in the numerical shock. To detect 
the shock we have used the basic check on the gradients of the numerical data. Define 

tj = max(\u(j ) - u(j - 1)|, | u(j + 1) - u(j)|) (5.1) 
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If tj > mox(3.0<j_2,3.0ij + 2,a), where a is chosen dynamically according to the structure 
of the shock wave, then we decide that the interval [xj_i, Xj+i] contains a discontinuity. 
Because of the global high order of the scheme, a falsely detected shock location in a smooth 
region of the solution would not destroy the overall accuracy of the scheme. 

To retain the spectral accuracy in the smooth region of the solutions, we apply high order 
exponential cut-off filters in (2.14). It is our experience that a very weak filter (i.e. high 
order) will suffice to get high accuracy in the smooth region. 

Time Discretization for Chebyshev Methods 


The time derivative in (4.1) is discretized with the Runge-Kutta method. We have used 
the third order Runge-Kutta method proposed in [17] which yields TVD (total variation 
vanishing) results if the spatial discretization is TVD. 

For periodical problems, the spatial derivative is approximated by Fourier trigonometric 
polynomials. When the solution is nonperodic, Chebyshev polynomials are used instead. A 
common difficulty, however, with Chebyshev methods is the stringent time step restriction. 
In general, the time step At has to be in the order of O(j^j-) where N is the order of the 
Chebyshev polynomials. This is a direct consequence of the clustering of Chebyshev collo- 
cation points near both boundaries. For many hyperbolic problems, this dense distribution 
of mesh points near boundaries is not necessary, as for most of the test problems in this 
paper. Recently in [12], a novel mesh transformation is proposed to relax the restriction 
of the Chebyshev methods on time steps. If x denotes the physical coordinate and £ the 
computational coordinate, the following transformation is considered in [12] : 


sin 1 ax 
sin -1 a 


M<i,l£l<i. 


(5.2) 


If & = cosi-^, 0 < i < N is the Chebyshev mesh in the space, then x, = £ sin(sin _1 a&) 
will be the corresponding mesh in the x— space. Because of the stretching nature of the 
transformation (5.2), mesh points x, will be more uniformly distributed in the physical x- 
space. A Chebyshev polynomial in the transformed £- space will be used to approximate the 
derivative with respect to £, and the derivative with respect to x will be computed as follows 


d d d£ a d 

dx d£ dx sin -1 acos(sin _1 a£) d£ 


(5.3) 


In our computations, we have observed an improvement of one magnitude in the time step 
with a = 0.999; at the same time, the resolution of the numerical method is also enhanced 
in the interior of the physical domain. We refer the reader to [12] for more details about the 
evaluation of the transformation. 
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Uniform Spectral Approximation for Discontinuous functions 

We consider the following piecewise C°° function 

{ — sin(2(x + 0.77 t)) + 1 — tt < x < — §7r, 

sin 2 x |*| < Itt, (5.4) 

J* " 1 - fx < X < 7T. 

</(x) is defined periodically on the outside of [ — 7r, 7 t] , thus g[x) has three discontinuities in 
the interval [0,27r]. The profile of g{x) is plotted as the solid line in Fig.2. Now, given the 
mesh value g(xi), x< = 0 < i < N, we use (2.20) and (2.21) to approximate g(x- + 1 ) 

and £g(xj), 0 < j < N respectively. In Fig. 1(a), (b), we plot the errors in function 
values and first derivative values respectively in the logarithm scale for different N’s, i.e. 
N = 32, 64, 128, 256. In all of the runs, we use a 16th order of exponential cut-off filter and 
the third order of ENO interpolation in (2.15). We clearly see the uniform convergence and 
spectral accuracy in the smooth parts of the function g(x). 

Linear Advection of Discontinuous Solution with Subcell Resolution 

In order to reduce the smearing in contact discontinuities by shock capturing schemes, 
Harten [10] suggested a sub-cell resolution technique to treat one- dimensional contact dis- 
continuities in the context of the cell-averaged ENO finite difference schemes. Later on, this 
idea was extended to the point-value version of the ENO finite difference scheme in [17]. We 
test the subcell resolution by our spectral methods. 

Consider the initial boundary value problem of the following linear hyperbolic equation 

Ut — U X} 

< u(x,0) = g(x), x € [0,2 tt] (5.5) 

u(0, t ) = u(2ir, t ) 

where ^(x) is defined in (5.4). 

In Fig. 2(a), (b), we plot the numerical solution for N = 200 after one and two period- 
icities, i.e. t = 2ir,t = Air. In both runs, we have used the 10th order exponential cut-off 
filters. 

Inviscid Burgers ’ Equation 

In this classic example of shock wave computation, we consider the initial value problem 
with sine wave initial conditions, 

' U t + (£ ), = 0 , 

< u(x,0) = a + /?sinx x 6 [0,2ir], (5.6) 

k u(0, t) = u( 2ir,t), 
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where a = 0.3, /3 = 0.7. 

The solution to (5.6) develops a shock discontinuity at time t = When a / 0, the 
solution consists of a moving shock wave after t = jj. As the exact solution for this problem 
can be obtained by iterative methods, we use this example to test the global accuracy of 
the scalar Algorithm I with m = 3 (i.e. the third order of ENO interpolation is used in 
(3.7)). In table 1, we list the global L\ error of the numerical solutions and the L\ error 
in the smooth region of the solution at time t = 2 and N = 32,64,128,256. In computing 
the global Li error, we exclude the occasional one transition point across the shock, and 
the smooth region is taken as the region which is 0.8 away from the shock location. The 
third column of Table 1 shows the global third order accuracy in L\ norm of Algorithm 
I, and the fifth column shows the increasing order of accuracy in the smooth region. In all 
of the runs, the time step has been chosen such that further decreasing of the time step 
does not improve the final accuracy. Therefore, the dominant errors come from the spatial 
discretization. In Fig. 3, solutions at time t — 2 and N = 32 are plotted (plus) against the 
exact 8olution(solid lines). In Fig. 4(a), the errors for N = 32, 64, 128, 256 at time t = 2 are 
plotted in logarithm scale. For a comparison, we plot the errors with the same parameters 
for the third order ENO finite difference scheme in Fig. 4(b). 


One- Dimensional Euler Equations for Gas Dynamics 

We solve the system of equations (4.1) with different initial data. 

1) Interaction of a shock wave and density waves 

We consider the following initial condition for (4.1), 

f pi = 3.857143, qi = 2.629367, Pi = 10.333333 -1 < x < -0.8, , 5 

( p T = 1 + esin57rx, q r = 0, P T = 1 —0.8 < x < 1, 

where e = 0.2. 

In Algorithm II, we choose the third order ENO interpolation and the 16th order 
exponential cut-off filters. In Fig. 5(a), we display the density profile for N = 200 at 
t = 0.36, the solid lines are the solutions obtained by the third order ENO finite difference 
methods with 800 points. For a comparison, we plot the solution obtained with the second 
order MUSCL schemes with N = 200 in Fig. 5(b). 


2) Sod’s problem and Lax’s problem 

We now consider the standard Riemann problem of (4.1) with the following initial data 
[16], 

a) Sod’s problem 
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b) Lax’s Problem 


(pi,q,,P,) = (1,0,1) -l<x<0, 

(p r ,q T ,P r ) = (0.125,0,0.10) 0 < x < 1; 


{ph qi, Pi) = (0.445, 0.698, 3.528) - 1 < x < 0, 

(/ Pr,q r >Pr ) = (0.5,0,0.571) 0 < X < 1. 

For the Lax’s problem, the sub-cell resolution has been applied on the linear degenerated 
field near contact discontinuities. In both of the problems, we use the third order ENO 
interpolation in Algorithm II and the 10th order exponential cut-off filters. In Fig. 6 
(a)-(c) the solutions of the Sod’s problems with N = 150 at time t = 0.4 are plotted, while 
Fig. 7(a)-(c) display the solutions of the Lax’s problem with N = 150 at time t = 0.26. In 
both cases, the solid lines are the exact solutions. 


Interaction of blast waves 


The initial data suggested in [20] to simulate the interactions of two blast waves is as 
follows 


u(x,0) 


” Uz, 

< U M 

. U H 


0 < x < 0.1, 
0.1 < x < 0.9, 
0.9 < x < 1, 


(5.8) 


where p L ~ p M = p R — l,qp = q\t = 9 r = 0,P L = 10 3 , Pm = 10 3 ,Pr = 10 3 . 

The solutions to this problem possess drastic fluctuations under the impact of inter- 
actions; it is a good test of the stability of Algorithm II. The complex structure of the 
solutions after the clash of two blast waves demands a stable high order method to capture 
the details of solutions. Unlike the finite difference methods, the spectral methods do not 


require exterior mesh points to treat boundary conditions. We apply characteristic boundary 
conditions on both boundaries. As the boundaries are treated as solid walls, we impose the 
condition that velocity variables vanish on both boundaries, i.e. qo = 0, q^ = 0. 

In Fig. 8(a)-(c) and Fig. 9 (a)-(c), we plot the solutions of the state variables with 
N = 300 at time t = 0.028 and t = 0.038, respectively. The former is an instant before the 
clash and the latter is one after the clash. The solid lines in both Fig. 8 and Fig. 9 are the 
solutions obtained with the third order ENO finite difference methods with 800 mesh points. 
In Fig. 10, the solutions of density with N = 400 are also plotted. 


Interaction between a Two-dimensional Shock Wave and a Rotating Vortex 


18 



The equations in consideration will be (1.1) with 


u = (p,m x , m y , E), (5.9) 

f(u) = q x u + (Q,P,0,q x P), (5.10) 

g(u) = q y u + (0, 0, P, q y P), (5.11) 

where q x , q v are velocity components in x- and y - directions respectively, m, = pq x and 
m y = pq y are x - and y - momentums respectively . P = (7 — 1)(I£ — \pq 2 ), q 2 = q x + q y - 


We apply one- dimensional Algorithm II on the fluxes f(u) and g(u) separately. The 
right and left eigenvectors for the Jacobian matrices can be found in [17]. 

The physical domain will be the rectangle [0,3] x [—1.5, 1.5]. A Mach 3 planar shock 
wave moves from the left to the right. A rotating vortex is initially located to the right side 
of the shock. As the time progresses, the shock will hit the vortex and interact with it. The 
shock front will be deformed by the interaction, and pressure waves are generated from the 
interactions. In the computations, we define the velocity fields of the vortex as those induced 
by two rotating concentric cylinders with radius r x and r 2 respectively, r x < r 2 . Initially the 
vortex is located at (x c ,y c ). The outside cylinder is stationary and the inside one rotates 
with the angular velocity u>. Let v(r ) be the radius velocity at a distance r from the center 
of the vortex, we have 

{ ur if 0 < r < r x , 

w r(^ + r ?) if 7 , i<r<r 2 , ( 5 - 12 ) 

0 if r > r 2 , 

where a = \ b = r\ - r\, r = J{x- x c ) 2 + ( y -y c ) 2 - We choose r x = 0.15, r 2 = 

r l r 2 V 

0.75, w = 7.5. 

Therefore the x - and y - velocities induced by this vortex at (x, y) will be 

?. = (5.13) 

r 

q y = i^£«(r), (5.14) 

r 

where x c = 2.25, y c = 0. 

The initial conditions for the simulation will be as follows: 


Pi 

= 3.857143, 


9x/ 

= 2.629367 if x < x 0 , 



= 0 , 

(5.15) 

Pi 

= 10.333333, 

(5.16) 
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and 


Pt — 1 > 

Qxr = Qx if X > X0> 

ftr = (517) 

P, = 1, 

where xo is the initial shock position, xo = 1.5. 

We impose characteristic boundary conditions on both left and right boundaries. A 
periodic boundary condition is considered in the y-direction and hence we are actually sim- 
ulating the interaction between an array of periodically distributed vortices and a planner 
shock wave. To relax the time restrictions of the Chebyshev approximation in the x - di- 
rection, we apply the mesh transformation (5.2) with a = 0.999. The shock has been made 
stationary by a translation in the mean flow direction. The second order ENO interpolation 
and the 10-th order exponential filter have been used in (4.10). 

Fig. ll(a),(b) are the contour plots of the pressure and density fields at time t = 0.4, 
while Fig. 12(a), (b) are the close-ups of the pressure and density at time t = 0.4. Fig. 13 is 
the pressure profile at time t = 0.4. 

Concluding Remarks 

Centered difference methods including spectral methods are efficient and accurate, while 
upwinding difference methods offer the advantages of sharp monotonic shock profiles. We 
have explored, in this work, the possibilities of blending the advantages of the ENO finite 
difference methods and the spectral methods. Numerical results have shown the robustness 
and feasibility of this approach at a small extra cost over the standard spectral methods. 
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Table 1: Global L\ error and L x error in the smooth region for the Burgers’ equations at 
time t = 2, the smooth region is defined to be 0.8 away from the shock location 


N 

Global 

Order 

Smooth region 

Order 

32 






1.49(-4) 


1.17(-4) 


64 


2.5 


4.3 


2.70(-5) 


5.86(-6) 


128 


2.9 


6.5 


3.70(-6) 


6.54(-8) 


256 


3.7 


10.0 


2.95(-7) 


6.36(-ll) 



23 

















Figure 2. Linear Advection of Discontinuous Solutions with Subcell Resolutions, N®200 
(a) time t » 2 7T ; (b) time t * 47 X 
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Figure 3. Solutions to the Inviscid Burgers' Equat 
Numerical Solutions (plus), Exact Soluti 
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Figure 4. Errors Co Che Inviscid Burgers' Equation for N = 32, 64, 128, 2 56 ; at j time c = 2 
(a) the Spectral Algorithm X, (b) the Third Order ENO Finite Difference Methods 





Figure 4. (Continued) 
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Figure 6. Solutions of the Sod's Problem with N = 150 at Time t = 0.4. (a) Density, (b) Velocity 
(c) Pressure. 





Figure 7. Solutions of the Lax’s Problem with N = 150 at time t ■ 0.26. (a) density, (b) velocity 
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Figure 7. (Continued) 
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Figure 8. (Continued) 
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Figure 8. (Continued) 









Figure 9. (Continued) 
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Figure 11. Interactions Betwwen a Two-Dimensional Planner Shock Wave and A Vortex, Shock Mach nurner 
(a) Density Contour with Level Value = 0.1, (b) Pressure Contour with Level Value = 0,2 
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Figure 11. (Co 
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Figure 12. Close-ups 
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Figure 12. (Conti 
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Figure 13. Same as in Figure 11, Density Solutions at Time 
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